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Abstract 

In this paper, we study the target tracking problem in wireless sensor networks (WSNs) using 
quantized sensor measurements under limited bandwidth availability. At each time step of tracking, the 
available bandwidth R needs to be distributed among the N sensors in the WSN for the next time step. 
The optimal solution for the bandwidth allocation problem can be obtained by using a combinatorial 
search which may become computationally prohibitive for large N and R. Therefore, we develop 
two new computationally efficient suboptimal bandwidth distribution algorithms which are based on 
convex relaxation and approximate dynamic programming (A-DP). We compare the mean squared error 
(MSE) and computational complexity performances of convex relaxation and A-DP with other existing 
suboptimal bandwidth distribution schemes based on generalized Breiman, Friedman, Olshen, and Stone 
(GBFOS) algorithm and greedy search. Simulation results show that, A-DP, convex optimization and 
GBFOS yield similar MSE performance, which is very close to that based on the optimal exhaustive 
search approach and they outperform greedy search and nearest neighbor based bandwidth allocation 
approaches significantly. Computationally, A-DP is more efficient than the bandwidth allocation schemes 
based on convex relaxation and GBFOS, especially for a large sensor network. 



I. Introduction 

A wireless sensor network (WSN) consists of a large number of spatially distributed sensors 
which are tiny, battery-powered devices, and have limited on-board energies. When properly 
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programmed and networked, WSNs perform different tasks that are useful in a wide range 
of applications such as battlefield surveillance, environment and health monitoring, and disaster 
relief operations. Dense deployment of sensors in the network introduces redundancy in coverage, 
so selecting a subset of sensors may still provide information with the desired quality. As shown 
in Fig. 1, the adaptive sensor management policies select a subset of active sensors to meet the 
application requirements in terms of quality of service while minimizing the use of resources. 
In this paper, we assume that the task of the WSN is to track a moving target in a given region 
of interest (ROI). Sensors receive observations from an object of interest and send quantized 
information to the fusion center over bandwidth limited channels. So the fusion center needs 
to distribute the available bandwidth among sensors using predictive information based on the 
target dynamics and the received sensor data. We consider a myopic (one-step ahead) scenario, 
where at a given time step, the fusion center only decides on the bandwidth distribution of the 
next time step. 



Fig. 1: System model for sensor and resource management based on feedback from recursive 
estimator. 

In the literature, there exist many sensor selection algorithms (see [1] and references therein). In 
[1], the sensor selection problem, an integer programming problem, has been relaxed and solved 
through convex optimization. One popular strategy for sensor selection is to use information 
driven methods [2], [3], where the main idea is to select the sensors that provide the most useful 
information, which is quantified by entropy or mutual information. The posterior Cramer-Rao 
lower bound (PCRLB) is also a very important tool because it provides a theoretical performance 
limit for a Bayesian estimator. As we have shown in our previous paper [4], for sensor selection, 
the complexity to compute the mutual information increases exponentially with the number of 
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sensors to be selected, whereas the computational complexity of Fisher information, which is 
the inverse of the PCRLB, increases linearly with the number of sensors to be selected. For 
target tracking in a bearing-only sensor network, a sensor selection approach which minimizes 
the PCRLB on the estimation error has been proposed in [5] and [6], where the selected sensors 
transmit either analog or quantized data to the fusion center. 

For the case where the fusion center receives quantized sensor measurements, given the total 

bandwidth constraint, R, at each time step during tracking, the fusion center should determine 

the optimal bandwidth distribution for the channels between the sensors and the fusion center 

which optimizes the target tracking performance in the WSN that consists of N sensors. This 

problem is more general than the sensor selection problem, because in the bandwidth allocation 

problem, the channel corresponding to each sensor could be assigned a different number of bits, 

while in sensor selection problems, a sensor is either activated or not to transmit its measurement 

under the constraint on the total number of selected sensors. The myopic bandwidth allocation 

problem can be solved by using an exhaustive search which enumerates all possible bandwidth 

distributions and decides on the solution that maximizes the determinant of the Fisher information 

matrix (FIM) which is the inverse of the PCRLB. Under Gaussian assumption, maximizing the 

determinant of the FIM is equivalent to minimizing the volume of the uncertainty ellipsoid [7]. 

/ R+N-l\ 

The search space of this problem is , which implies that explicit enumeration 

V J 

of all the solutions is computationally prohibitive for large and R. Therefore, computationally 
efficient suboptimal methods are required. In [8], the generalized Breiman, Friedman, Olshen, 
and Stone (GBFOS) algorithm has been employed for dynamic bandwidth distribution for target 
tracking which significantly outperforms a static equal bit allocation scheme in terms of tracking 
performance. But still, as we show later in the paper, the GBFOS algorithm may become 
computationally costly with increasing values of A^. 

Dynamic programming (DP) [9] solves the resource allocation problems by breaking them 
down into simpler steps. For a scalar- valued parameter estimation problem, a DP recursion can 
be easily formulated to find the optimal bandwidth distribution at each time step by maximizing 
the Fisher information due to the fact that the total Fisher information is the summation of each 
sensor's individual Fisher information. For target tracking, even though the Fisher information 
is in a matrix form and the objective is to maximize the determinant of the FIM, we can still 
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formulate a DP recursion which would yield a suboptimal solution. We refer to this scheme as 
approximate DP (A-DP), which is computationally very efficient since its complexity increases 
linearly with N. 

In our preliminary work [10], we compared the performances of dynamic bandwidth allocation 
approaches based on A-DP, GBFOS and greedy search. Motivated by the sensor selection 
method presented in [1], in this paper, we first formulate the bandwidth allocation problem as a 
constrained optimization problem with binary-valued decision variables and equality constraints. 
We then relax and solve the problem optimally using Newton's method by replacing the Boolean 
variable, m G {0, 1}, which represents whether or not the quantized measurement of sensor 
i is transmitted to the fusion center in m bits, with its convex counterpart gj „i e [0, 1]. Using 
the idea of probabilistic transmission for bandwidth management [1], [11], we treat gj „i as the 
transmission probability, that is at a given time, sensor i transmits its decision to the fusion 
center in m bits with probability ,„ G [0, 1]. Therefore, the convex relaxation based bandwidth 
allocation method meets the bandwidth constraint in an average sense and introduces a weak 
constraint on bandwidth availability. We compare the bandwidth allocation schemes based on 
convex relaxation, A-DP, GBFOS and greedy search in terms of their mean squared error and 
computational load under different process noise parameters. Simulation results show that convex 
relaxation, A-DP and GBFOS yield similar tracking performance, which is also similar to that 
of the optimal bandwidth allocation scheme based on exhaustive search. Among these three 
suboptimal schemes, A-DP has the least computational load, when the sensor network is large. 

The rest of the paper is organized as follows. In Section II, we introduce the target tracking 
problem, and describe the optimization of the quantization thresholds and particle filtering in 
target tracking. In Section III, we describe the bandwidth distribution schemes based on convex 
relaxation, A-DP, GBFOS and greedy search. In Section IV, we present numerical examples and 
compare the performances of the considered bandwidth distribution schemes in terms of their 
computational load and MSEs. Finally, we conclude our work in Section V and discuss some 
future research directions. 

II. Target Tracking in Wireless Sensor Networks 

The problem we seek to solve is to track a moving target using a WSN where N sensors 
are grid deployed in a square surveillance area of size 6^. The assumption of grid layout is not 
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necessary but has been made here for convenience. Target tracking based on sensor readings can 
be performed for an arbitrary network layout if sensor placements are known in advance. All the 
sensors that are assigned bandwidth report to a central fusion center, which estimates the target 
state, i.e., the position and the velocity of the target based on quantized sensor measurements. 
We assume that the target (e.g., an acoustic or an electromagnetic source) emits a signal from the 
location (xt,yt) at time t. We assume that the target is based on flat ground and all the sensors 
and target have the same height so that a 2-D model is sufficient to formulate the problem. 

At time t, the target dynamics are defined by a 4-dimensional state vector xt = 
[xt Ut it VtY where Xt and iji are the target velocities in the horizontal and the vertical 
directions respectively. Target motion is defined by the following white noise acceleration model: 



(1) 



where F models the state dynamics and Vt is the process noise which is assumed to be white, 
zero-mean and Gaussian with the following covariance matrix Q. 
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(2) 



In (2), V and p denote the time interval between adjacent sensor measurements and the process 
noise parameter, respectively. It is assumed that the fusion center has perfect information about 
the target state-space model (1) as well as the process noise statistics (2). 

The target is assumed to be an acoustic or an electromagnetic source that follows the power 
attenuation model provided below [12]. At any given time t, the signal power received at the 
sensor i is given as 



''i,t 



adit 



(3) 



By adopting this model, we prevent the receiver amplifier from saturation and the regularity 
conditions for PCRLB hold when the target is very close to a sensor. In Eq. (3), Pq denotes the 
signal power of the target, n is the signal decay exponent and a is a scaling parameter, dj ^ is the 
distance between the target and the sensor, di^t = a/ {xi — XtY + {jji — UtY, where {xi,yi) 
are the coordinates of the i^^ sensor. Without loss of generality, a and n are assumed to be unity 
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and 2, respectively. At time t, the received signal at sensor i is given by 

Zi,t = ai^t + rii^t (4) 

where t is the noise term modeled as additive white Gaussian noise (AWGN), i.e., rzj t ~ 
A/'(0, 0"^), which represents the cumulative effects of sensor background noise and the modeling 
error of signal parameters. 

Rather than transmitting analog sensor observations to the fusion center, transmitting a 
quantized version of sensor measurements decreases the amount of communication and therefore 
reduces the energy consumption. A sensor measurement Zi^t at sensor i is locally quantized before 
its transmission to the fusion center using Ri t bits. Let = [Ri^t, ■ ■ ■ , Rw.t] be the vector 
of quantization rates used by the N sensors in the network. For the bandwidth distribution 
problem, at each time step of tracking, Ri t can take values m, m G {0, 1, . . . , i?} where R 
is the maximum number of bits to be transmitted to the fusion center collectively by all the 
sensors. Let Lm 4 2™ — 1 be the number of decision intervals for transmitting m bits to the 
fusion center and Di t be the m-bit observation of sensor i quantized with rate i?j t = m at time 
step t, then 

-oo < Zi^t < VT 

1 vT < Zi,t < vT 



D 



•Lt 



(5) 



^ L„, - 1 ?7lL_i < Zi^t < oo 

where rf = [77^ 77™ ... r]^^] with 77™ = — 00 and r/^'^^ = 00. The quantization thresholds are 
assumed to be identical at each sensor for simplicity. We explain the selection of the quantization 
thresholds for each data rate R^ t = m later in this section. Given Xj and m, it is easy to show 
that the probability of a particular quantization output / is, 

P(A, = /|x„ R., = m) = Q (^^^^^) - Q (6) 

where Q{.) is the complementary distribution function of the standard Gaussian distribution with 
zero mean and unit variance, 

Q{x)= j ^expf-Mrft (7) 
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At time t, let the fusion center receive the data vector = [-Di,t, . . . , DN,t] from the sensors 
with the corresponding quantization rate vector Rt = [Ri,t, • • • , RN,t], then 

N 

p(Dt|xi, Rt) = JJp(A,t|xt, Ri,t) (8) 
1=1 

where we assume p(A,t|xf, = 0) = 1. 



A. PCRLB with quantized data 

Let p(Dt, X() be the joint probability density of and x^, and be an estimate of x^ at time 
step t. Based on the received data Dt quantized with rate vector R^, and the prior probability 
distribution function of x^, p(xt), the PCRLB on the mean squared estimation error has the form, 

E{[±t- Xi][xi - ^tf\Kt} > Jr'(Rt) (9) 

where Jt(Rf) is the 4x4 Fisher information matrix (FIM) with the elements 

logp(Dt,xt|Rt 



z,jG{l,...,4} (10) 



dxt{i)dyit{j) 

where Jt (Rt) (i, j) denotes the i*'^ row, j*'* column element of the matrix Jt(Rt) and Xt(z) denotes 
the i*'* element of vector Xf. Let V^J = Vxt V^^ denote the second order partial derivative operator 
with respect to x^. Using this notation, (10) can be rewritten in a more compact fashion as, 

Jt(Rt) = E [-V^:iogp(Dt,xt|Rt)] (11) 

Since p(Dt,xt|Rt) = p(Dt |xt, Rt)p(xt), Jt(R.t) can be decomposed into two parts as, 

Jt(Rt)=jf(Rt) + jf (12) 

where 

jf(Rt) ^i?p(D,|x,Mx.) [-V^:iogp(D,|x„Rt)] 
Jf ^E,(xO [-V^qogp(xt)] 

Note that Jf(Rt) represents the Fisher information obtained from the data averaged over the 
prior distribution p(xt) and Jf represents the a priori Fisher information. -E'p(Dt|xt)p(xt)[-] and 
£'p(x()[.] denote expectations with respect to p(Dt|xt)p(xt) and p(xt) respectively. 
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Given the vector of quantization rates Rt = [Ri,u ■ ■ ■ ,RN,t] and using (8) in (12), the data 
part of the Fisher information can be written as, 



Jf{Ri,u...,RN,t) =[ ^p(D,|x„R,) [-^x:iogP(Di|xi,Ri)]p(xi)t/xi (13) 

i = l >^Xt J 

For a given Xj, let us define Jfj(i?j tlxj, as the Fisher information of sensor i, 

(14) 



= Yl {-v^:iogp(A,t = /|xt,i?,,t)p(A,i = /|xt,i?,,t)} 



1=0 



Then combining (13) and (14), sensor i's contribution to the Fisher information Jf^^(i?j t) can 
be stated as. 



(15) 



Given Rj, the Fisher information at time t can be written as, 

N 



J,(R,) = J]jg(i?,,) + J 



(16) 



i=l 



From (14), after straight-forward calculations, the (1, 1) term of Jf^{Ri^t\^t) can be derived as. 



E 



d'^ logp(A,f|xt 
dx'i 



■i,t I 



p{Di^t = l\^t,Ri,t) \ 



dxt 



(17) 



The rest of the terms can be derived similarly. Using the procedures similar to [13], Jf^{Ri^t\^t) 
can be obtained as follows. 



Jit{Ri,t = m\yit) = n Ki^t{m,Xi,yi,Xt,yt)-r , ,„ , „ 

(1 + ad 



(18) 



(xi-xt)^ {xi - xt){yi - yt) 

{xi - xt){jji - yt) ijjt-ytf 
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where 



— 2^ — 



2 




1 



e 



— e 



> 



(19) 



V 



Detailed derivation of (18) can be found in the Appendix. Note that in (18) and (19), da and 
Oj t are functions of the sensor location {xi,yi) and target location (xt,yt)- 

B. Optimization of Quantization Thresholds 

The Fisher information and hence the PCRLB are functions of the quantization thresholds 
corresponding to each data rate Ri^t = rn. Thus, the quantization thresholds should be designed 
to achieve better estimation accuracy. An algorithm to obtain the optimal quantization thresholds 
that minimizes the variance of the estimation errors has been proposed in [13]. If we assume 
that {xi,yi) and {xt,yt) are uniformly distributed in a region, we can minimize the sum of 
two diagonal elements of the CRLB matrix, after averaging the CRLB matrix over all the 
random parameters which may result in a large computational load since it requires a multiple 
fold integration. To alleviate this problem, some alternative methods to design the quantization 
thresholds were developed in [13]. 

Note that all the information about [xt,ytf^ is contained in sensors' signal amplitudes (aj t)'s. 
If all the signal amplitudes can be recovered from their quantized data Di^t accurately, an accurate 
estimate of [xt, yt]'^ can be obtained. In this paper, we use the Fisher information based heuristic 
quantization method [13] which maximizes the Fisher information about the signal amplitude 
tti^t contained in the quantized data Di^f We define Fa{r}\xi,yi, Xt,yt, Ri,t = rn) as the Fisher 
information of the signal amplitude contained in quantized m-bit data, Di t, using a threshold 
77. Note that a^ ^ is a function of di t for fixed Pq, a and n as defined in (3). Then given 
Ri,t = rn, sensor location {xi,yi) and source location {xt,yt), it has been derived in [13] that 
Faivl^i^Ui^ ^t,yt, Ri,t = m) = 4:Ki t{m,Xi,yi,Xt,yt)- The Fisher information based heuristic 
quantization method [13] finds the decision thresholds that maximize 



FairilRi^t = m) = E[ 



Va^;* log p{Di^t\ai^t{xi, y^, Xt, yt))] 



(20) 
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where u 



d}f and the Fisher information about the signal amplitude is averaged over 



the probability density function of u, p{u), under the assumption that {xi,yi) and {xt,yt) 
are independent and identically distributed and follow a uniform distribution U[—b/2,b/2]. 
Derivation of p{u) and other details of this quantizer design approach can be found in [13]. 
We assume that the decision thresholds of each quantization rate are identical at each sensor. 
The quantization thresholds of each possible quantization rate is optimized offline and can be 
stored at each sensor before the WSN operation. 

C. Particle Filtering with Quantized Data 

It is known that Kalman Filter provides the optimal solution to the Bayesian sequential 
estimation problem for linear and Gaussian systems. In nonlinear systems, the extended Kalman 
filter (EKF) can be used to provide a suboptimal solution by linearizing the nonlinear state 
dynamics and/or nonlinear measurement equations locally. However, it has been shown [14] 
that, even for linear and Gaussian systems, when the sensor measurements are quantized, the 
EKF fails to provide an acceptable performance especially when the number of quantization 
levels is small. Therefore, we propose to employ a particle filter to solve the Bayesian sequential 
estimation problem. 

Let Dit = [Di,...,D(] be the received sensor data up to time t which are obtained 
according to the data rates Ri;^ = [Ri, . . . , Rf]. In particle filtering, the main idea is to find 
a discrete representation of the posterior distribution p{xt\T)i;t) by using a set of particles 
{x^; s = 1, . . . , Ng} with associated weights {wf; s = 1, . . . , A^^}- The posterior density at 
t can be approximated as. 



where Ng denotes the total number of particles. In this paper, we employ sequential importance 
resampling (SIR) particle filtering algorithm [15] to solve the nonlinear Bayesian filtering 
problem. In Algorithm 1, we provide a summary of the SIR based particle filtering rather than 
discussing the details. Note that Tg in Algorithm 1 denotes the number of time steps over which 
the target is tracked. A more detailed treatment of particle filtering can be found in a wide 
variety of publications such as [15]. 




(21) 



s=l 



October 26, 201 1 



DRAFT 



JOURNAL OF MbX CLASS FILES, VOL. X, NO. X, SEPTEMBER 20XX 



11 



Algorithm 1 SIR based Particle Filtering for Target Tracking 

Set t = 0. Generate initial particles Xq ~ p(xo) with ,Wq = N~^. 
while t < Ts do 

(Al.l) = Fx.1 + Vt (Propagating particles) 
(A1.2) |)(x,+i|Di,) = ijEf=i5(x*+i -x?^i) 

(A1.3) Bandwidth Allocation: Decide Ht+i and obtain sensor data Dj+i 
(A1.4) wf+i oc p(Dt+i|xf+i,Rt+i) (Updating weights) 
w^_^^ = (Normalizing weights) 

^ v^A^s s s 

— 2^s=i ^t+i^t+i 
(A1.5) {x,%i,iV7i} =Resampling(x,%i,u;f+i) 

(A1.6) t = t + l 

end while 



In Algorithm 1, j>(Dt+i|x^^]^, R^+i) is obtained according to (6) and (8). Resampling step 
avoids the situation that all but one of the importance weights are close to zero [15]. 

By using equations (10) to (16), at time t, one can compute the PCRLB on the estimation error 
and the corresponding FIM, for a given bandwidth allocation scheme Ht and prior distribution 
p{xt). For the bandwidth allocation problem, at time t, from (A1.2), we first generate the prior 
p{xt^i\Y)i,t) using data received up to time t. 

Under the Gaussian assumption, maximizing the determinant of the FIM is equivalent to 
minimizing the volume of the uncertainty ellipsoid [7]. Therefore, we determine bandwidth 
allocation scheme for time t + 1, Rt+i, by maximizing the determinant of the Fisher information 
about xt_|_i as, 

max/j,^,+,,...,R^^^^^ det(Jt+i(Rt+i)) (22) 

N 

S.t. ^ Ri,t+i = R 

i=l 

The fusion center then informs the sensors about Rt+i and sensors transmit their quantized 
measurements D^+i accordingly. The Fisher information, Jt+i(Ri+i) is written as 

Jt+l(Rt+l) = ^p(D,+i,x,+i|Di:i,R,+i) [-V^;+Mog|)(Dt+i,Xt+i|Di:(,Rt+i))] 
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Following the derivations from (10) to (16), the Fisher information, Ji+i(Rt+i), is obtained as, 

N 

Jt+i(Rt+i) = ^ J(^i(i?i,t+i) + Jt+1 (23) 

i=l 

Using the particle approximation, 

. Ns 

p(x,+i|Di,) ^ ^5(x,+i - x,Vi) (24) 



s=l 



Jf^]^(i?j^t+i) is found from. 



J^i(-Ri,t+i) = Jf+i(^i,i+i|xt+i) (25) 

s=l 

As in (12), 3f_^^ = £'p(xt+i|Di:t)[— Vxt+J logp(xt+i|Di:j)] has been defined as the prior Fisher 
information of x^+i. According to (24), p{xt^i\Di,t) has a non-parametric representation by a 
set of random particles with associated weights, so it is very difficult to calculate the exact Jf^^ 
[16]. Instead, we use a Gaussian approximation such that j>(xf+i|Di:t) ^ J\f{fi^^-^^, St+i), where 

1 



s=l 



and 



1 



^ s=l 



Given the Gaussian approximation, it is easy to show that J^^^ = ^t+i- 



III. Dynamic Bandwidth Allocation for Target Tracking 

An exhaustive search can be employed to find the optimal bandwidth distribution which 

maximizes (22). For a network of N sensors and bandwidth constraint R, there are a total of 

R + N — 1 \ (jv+R-i)' 

— (jv-ivm possible bandwidth distribution solutions. For large N and R, such 

N-1 J ^ 

an exhaustive search may not be feasible in real time. Therefore suboptimal but computationally 
more efficient algorithms are required which we explore in this section. 
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A. Convex Optimization Based Dynamic Bandwidth Allocation 

In this paper, we use the log determinant of the FIM as the objective function for resource 
management. Using Boolean variables ^ G {0, 1}, the bandwidth allocation problem can be 
explicitly formulated as follows, 

R N 



\m=0 i=l 



maxq^^^ logdet(Jt+i(q^+i)) = logdet | ^ ^ qi^rnJ^+ii^^Lt+i = m) + J^^ j (26) 

R 

subject to ^ qi^rn = 1 i G {1, . . . , A^} 



m=0 
R N 



m=0 i=l 

gi,„G{0,l} mG {0,l,...,i?} iG{l,...,iV} 

In the above formulation, q^^^ = [gi,o, ^2,0, • • • , g7v,o, • • • , ^2,^, • • • , g^v,/?]'^ denotes the 
bandwidth allocation scheme for time t + 1 where qi^m = 1 when sensor i transmits its 
measurement in m bits and Jf^_^_^{m) is the corresponding FIM of sensor i. Note that we drop 
the time index t + 1 from the elements of vector q^+i to simplify the notation. All constraints 
are equality constraints where the first constraints guarantee that each sensor can transmit 
using only one of the quantization rates. If m = is selected, the quantized measurement of the 
sensor is not transmitted to the fusion center. The (N + 1)*^ constraint ensures that the sum rate 
constraint is satisfied and the last N{R + 1) constraints restrict „ to be Boolean. 

Similar to the convex relaxation approach presented in [1], the last N(R+ 1) constraints can 
be relaxed by replacing the Boolean variables ^ G {0, 1} with their continuous counterparts, 
(!i,m G [0, 1]. Then the problem becomes 

/ R N \ 
max^^^^ logdet(Jt+i(qt+i)) = logdet ^ ^ gi,^jf^4^i(i?i,i+i = m) + J^^ (27) 



\in=0 i=l 



R 

su'j 

m=0 
R N 



bject to ^ qi^m = 1 i e {1, . . . , N} 



^ ^ m qi^rn = R 

m=0 i=l 

0<q^,„^<l mG {0,l,...,i?} zG{l,...,iV} 
We can further relax the problem (27), by removing the last N{R + 1) constraints and including 
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them into the objective function. Then, the new cost function to be minimized becomes as, 



log det I ^ ^ qi^rnJ 



(28) 



R N 



D 



m 



Jt+l I + ^ X] ( log i%rn) + log (1 - qi,m) j > 

J m=0 1=1 ^ ^ J 



\m=0 i=l 

where 0(qt+i) is a convex function of the decision variables gj „ [17]. The additional summation 
term in the objective function forces m to be in the interval [0, 1]. r is a positive parameter 
that controls the quality of the approximation. 
Let us define. 



/l 
1 



A 



1 
1 



... 1 
... 



\ 





\ 



1 
1 1 



1 ... 
1 ... R 



1 
R R J 



qt+i 

and 



Qi,o Q2,o ■ ■ ■ Qn,o Q2,i 



1 1 ... 1 R 



Qn,i 



<1n-i,r Qn,r 



Then, the first + 1 equality constraints of (27) can be represented in a matrix form as, 

Aqt+i = b 

Finally, we have the following convex optimization problem, 

minq,+, 0(qt+i) 
subject to Aqt+i = b 



(29) 



which can be solved efficiently and optimally using Newton's method. The underdetermined 
system Ac^t+i = b has infinite number of solutions but there is only a subset of solutions 
which are feasible satisfying < q^+i < 1 where and 1 are the all zero and all one vectors 
respectively. The Newton method starts with a feasible solution, so we formulate the following 
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linear optimization sub-problem to find an initial feasible solution, 

R N 



min^^^j ~ 5Z 5Z ^''"^ (^^^ 

m=0 i=l 

subject to Ac{t+i = b 

0<gv;„<l me{0,...,R} ie{l,...,N} 

The optimality conditions for (29), which is named as the Karush Kuhn Tucker (KKT) system, 
is written as [17], 



(31) 



A 

where AZ is the Newton Step, u is the the optimal dual variable, V q^^^cf) is the gradient vector, 
and V^'^|0 is the Hessian matrix of with respect to the decision vector qt+i. 

In order to solve the system given in (31), first we need to compute the gradient vector, 
and the Hessian matrix, V|+|0, with sizes N{R + 1) x 1 and N{R + 1) x N{R + 1) 
respectively. Let us start by computing the gradient. First, the (z, m)*^ element of the gradient 
vector is, 

(Vg,^,0),,^ = (32) 

d \ f R N \ Q R N 

—— I - log det ^ gi,„jg+i(m) + J^^ - r—— ^ ^ (log (g^,^) + log (1 - g^,™)) 

I \m=0 i=l / m=0 

Let X be an invertible matrix and x be a scalar. Using the property 

aiogdet(X) _^ r^„i9X 



dx I dx 

and the definition 

(R N 
EE^^--JS+iM + jf+i 
m=0 i=l 

each element of the gradient vector is obtained as 

{V^.^A\ra= -tr{W-ijg+i(m)}-^ + — ^ (33) 
In order to compute V^^^|0, for E {1,2,..., A^} and m,m* E {0, 1, ... , R}, we define 
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Using the properties 



dx \ dx 



and 

dx dx 

we get 

^i^^,m^U^,^) = tr{W-^J,^,+iK)W-ijg+i(m)} 
Finally the Hessian matrix is obtained as 

(V|::0) =^ + r drag + 7^^; • • • 5 J- + TT^^) (35) 

Having obtained Vg^^^^, V^*^|0 and A, the block elimination method [17] can be used 
to solve the KKT system, which is summarized in Algorithm 2. (For the details of block 
elimination method see Section 10.4. in [17]). In order to compute the complexity of the 

Algorithm 2 Solving KKT system by block elimination 
(A2.1) Form (V|+»-U^ and (V|^»-i(V^,^,0). 
(A2.2) Form S = -A{Vf+l(j)y'A^ 

(A2.3) Determine uj by solving Scu = ^(V|+|</))-^(Vg,+i</)). 
(A2.4) Determine AZ by solving (V|+|0)"^AZ = A^oo - 



convex relaxation based bandwidth allocation method, we ignore the complexity for computing 
(Vgt_^i0) and (V?'^|0). The cost of block elimination to solve (29) is dominated by the Cholesky 
decomposition of V^|^|0 which is used to find (V?*^|</))^^. Let us define V = N(R + 1). Then, 
in order to compute the Cholesky decomposition of the matrix (V^|^|0), we require a total of 
|(V^ — V) summations and multiplications and |(3V^ — 3V) divisions [18]. Thus the complexity 
of bandwidth allocation based on convex optimization increases with 0{N^{R + 1)"^). 

At each iteration of Newton's method, the solution vector q^+i is updated by q^+i = q^+i + 
sAZ where s E (0, 1] is the step size obtained by the backtracking line search method [17]. We 
stop the Newton iterations when the Newton decrement, A = (— Vg^^^^-^AZ)^/^, is less than 
some predefined precision. The summary of the Newton's method is presented in Algorithm 3. 
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Algorithm 3 Newton's method for the bandwidth allocation problem 
Find a feasible starting point q^+i from (30) and set precision e > 

Repeat 

(A3.1) Compute the Newton Step AZ from (31) and Newton decrement A 

(-V,,,,0^AZ)V2 

(A3. 2) Choose Step size s by backtracking line search. 
(A3.3) Update q^+i = q^+i + sAZ 

(A3.4) Quit if \'^/2 < e, then q*^i = q^+i, else go to Step (A3.1). 



The Probabilistic Transmission Scheme: From the optimal solution of the relaxed problem 
^l^i, the bandwidth distribution for the next time step needs to be determined. For the sensor 
selection problem, the authors in [1] employed a simple scheme in which k sensors are selected 
out of N sensors by first sorting cil_^_i in descending order and then setting k largest elements of 
ql^i to one. For the bandwidth distribution problem, a similar solution is to sort the probabilities 
Qi^rn in descending order and then to assign 1 starting from the largest probability until the 
bandwidth constraint is satisfied. However, in this paper, we consider a randomized scheme 
similar to the ones used in [1] and [11]. Since the elements of qj'^^ are within the range [0, 1], 
we can consider each „ as the transmission probability of sensor i, transmitting information 
in m bits. Instead of putting a strict bandwidth constraint, i.e. J2f=i Ri,t+i = the probabilistic 
transmission puts a weak constraint on the bandwidth availability and ensures that the sensors 
on the average transmit R bits to the fusion center. We present a numerical example on the 
probabilistic transmission scheme in Section IV. 



B. Approximate Dynamic Programming based Bandwidth Distribution 

In this section, we present the bandwidth allocation algorithm based on A-DP which will be 
shown to provide near optimal solution but require much less computation time than the convex 
relaxation approach. Note that the Fisher information matrix can be expressed as the summation 
of each sensor's individual Fisher information matrices as defined in (23). In this section, we 
formulate an approximate DP recursion in tracking applications where we can maximize the 
Fisher information by maximizing its determinant subject to the bandwidth constraint. 
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Typically DP involves progression along time. But in our problem formulation the DP 
progresses across sensors and is executed at each time step of tracking to determine the bandwidth 
allocations of the next time step. Since A-DP is performed at each time step, for simplicity, the 
time index t+1 for Fisher information matrix is dropped. Instead an index for the stages in DP 
is adopted. Let Jjv = Jt+i and Aj(i?j) = J^^(i?j) be the reward in terms of Fisher information 
when sensor i quantizes its measurement in Ri bits (Ri E {0,1,..., R}). While constructing the 
DP trellis, the bandwidth distribution problem is first divided into N+1 stages which correspond 
to sensors and a termination stage. We define the state of a stage as the remaining bandwidth 
for the usage of sensor i. So each stage has R + 1 states associated with it. The bandwidth 
allocation chosen at any sensor (stage) determines the feasible states at the next sensor. An 
example DP trellis is shown in Fig. 2 with A^ = 6 and R = 3 which implies a total of 7 stages 
and 4 states in the DP trellis. As an example, sensor 1 is at state r = 1 means 2 bits have already 
been used by the other A^ — 1 sensors and 1 bit is available for sensor 1. Then, sensor 1 can 
only take the action Ai(l) and the DP goes to the termination stage (stage 0) which has only 
the bit available state. 

For such a DP trellis, we have, 

Jjv = Ajv(i?jv) + {An-i{Rn-i) + ■■■ + Ai(i?i) + Jo} 

Ji =Ai(i?i) + Jo (36) 
where Jq = and = According to the matrix determinant lemma [19], 

det(X + A) = det(X + AI) = det(X) det(I + X^^A) 
With X = Jj_i, A = Aj(i?j), and I being the identity matrix, we have 

log I det(J7v)| = log I det(Jjv-i)| + log | det [l + J^^_iA;v(i2^)] | 

log I det(Ji)| = log I det(Jo)| + log | det(I + Jq ^Ai(i?i))| 
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We can maximize det(JAr), by maximizing log | det(JAr)|. The DP recursion at each stage 
is formulated as follows: the trellis starts from Jo(0) and for the first stage (i = 1) and for allr G 
{0,1,. 



log 



det[Ji(r)] 



log 



det[I + Jo-'(0)Ai(r)] 



+ log 



det[Jo(0)] 



Then for all the intermediate stages i E {2, . . . , N — 1} and for all r G {0, 1, . . . , R}, 



log 



det[J,(r)] 



(37) 



(38) 



max < log 

fc=0,l,...,r 



det[I + Jr\(r- A;)A,(/c)] 



log 



det[Ji_i(r - k)] 



Finally for the last stage i = N, 



log 



det[J7v(i?)] 



(39) 



max 

k=OA,....B. 



log 



+ log 



det[JN-i{R- k)] 



In (37), (38), and (39), the reward of sensor z's transmission in Ri bits depends not only on 
Aj(i?j) but also on the FIM of the previous stage Jjl\. So at each stage i, the FIM, Ji(r), which 
has the maximum determinant should be stored in a memory for its use at the next recursion. 
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Note that the proposed A-DP may not yield the maximum matrix determinant at the final stage. 
The suboptimality of the A-DP recursions is discussed later in this section. 

We analyze the computational complexity of A-DP in terms of number of matrix summations. 
Note that the number of element-wise summation is a scaled version of number of matrix 
summations. The first stage needs R matrix summations to compute the FIM at all states. For 
all the intermediate stages, at each state r, (r G {1, . . . , i?}), r different matrix summations are 
required to find the FIM with the maximum determinant. Finally at stage N, A-DP again needs 
over R matrix summations in order to maximize the determinant of det Jat. So the A-DP totally 
searches over, 

R 



R+{N -2) 



r 

r=l 



+ R = 2R+{N-2f^^^^^ 



2 

matrix summations which is linear in and quadratic in R. 

Suboptimality of the A-DP recursions: For a given state of a stage, we choose the path with 
the maximum determinant of the FIM and dismiss all the other paths arriving at this state. The 
proposed DP recursions would yield the optimal solution to maximize the determinant of the 
FIM, if the following property were satisfied, 

if det{j'} > det{j"} (40) 
then det{A + j'} > det{A + j"} 

for some positive semidefinite matrices j', j" and A. Unfortunately, the above property is not 

/ 1 \ „ ( I -0.1 

necessarily true. Consider the simple example, J = and J = 

\^ 1 y \^ -0.1 1 

( I 

where det{J } > det{J }. Let A = . Then det{A + J } < det{A + J }. 

At each stage of the DP, we only store the FIM with the maximum determinant. Therefore, 
the final solution obtained by the DP recursions becomes suboptimal since not all the feasible 
solutions are enumerated. 



C. Existing Suboptimal Bandwidth Distribution Methods 

In this section, we review some existing suboptimal methods that are suitable for solving the 
bandwidth allocation problem in target tracking applications. 
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1 ) GBFOS Algorithm: This algorithm has been first proposed in [8] for dynamic bandwidth 
allocation in target tracking. The GBFOS algorithm starts by assigning the maximum number 
of bits, R to each sensor in the network and then reduces the number of bits one bit at a time 
until the sum rate constraint is satisfied. The GBFOS algorithm can be stated as in Algorithm 
4. Note that in order to simplify the notation, we drop the time index t + 1 in the algorithm. 
As shown in Step (A4.1), at each iteration, the GBFOS algorithm searches the sensors and 
reduces the bits of the sensor by one which ensures the minimum reduction of the determinant 
of the FIM. An efficient implementation of the GBFOS algorithm and its complexity analysis 
can be given as follows: 

Let us define J* as the FIM after the z*^ iteration, and Ak{Rk) — Jt+i{Rk) as sensor A;'s 
contribution to the FIM using R^ bits. In the beginning, we need to calculate J° = + 
Ai{R) + A2{R) + .... + An{R). As a result, totally N matrix summations are needed. At the 
i-th iteration, where i = 1, . . . , {N — l)R, there are at most different ways to reduce 1 bit. 
Assuming one particular solution is to reduce one bit at the A;-th sensor, {k E {1,2,..., A^}), 
and J*(/c) = J*^^ + Ak{Rk ~ 1) — Ak{Rk), which requires two matrix summations. Hence at 
each iteration, at most 2N matrix summations are required. At the end of the i*'* iteration, we 
store J^{k) with the maximum determinant. In summary, we need at most A^ + 2N{N — 1)R 
matrix summations, which is quadratic in A^ and linear in R. Note that this is an upper bound 
on complexity. 

2) Greedy Algorithm: Basically, greedy search is the reverse of the GBFOS method which 
makes the algorithm much faster. The greedy algorithm can be stated as in Algorithm 5. The 
greedy algorithm starts by assigning bits to each sensor in the network and then increases the 
number of bits one bit at a time until the sum rate constraint is satisfied in R iterations. At each 
iteration, greedy algorithm searches the A^ sensors and a single bit is added to the sensor which 
maximizes the determinant of the resulting FIM. 

The implementation of greedy search and its complexity can be stated as follows: At the first 
iteration, there are A^ different ways to add 1 bit. For the fc-th way of adding 1 bit, J^(fc) = Jt+i + 
Afc(l). Then we set = max^ det(J^(A;)). Hence, totally A^ matrix summations are required at 
the first iteration. At the i-th iteration, for i = 2, . . . , R, there are still A^ different ways to add 1 
bit. For the A;-th way of adding 1 bit, J'{k) = J*~^ + Ak{Rk + 1) - Ak{Rk)- Note that Ak{Rk) 
could be a zero matrix since Rk could be zero. Therefore, for each iteration, at most a total of 
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Algorithm 4 GBFOS - Bandwidth Distribution Algorithm 



Set Ro = [Ri = R,...,Rn = R] and J° = Jf^^ + Ai{R) + A2(i?) + .... + Ajv(i?). 
FOR t = 1 : {N - 1)R 
(A4.1) FOR k=l: N 
IF Rk>0 

Reduce one bit from sensor k and compute det(J*(A;)) where 
ji(k) = + AkiRk - 1) - AkiRk). 
ENDIF 
ENDFOR 

(A4.2) V k with Rk > 0, find the sensor p* for which det{3\k)) is the maximum: 

p* = arg max det(J*(A;)). 

k where Rk>0 

(A4.3) Decrement Rp* = Rp* — 1, update Rj = [Ri, . . . , Rp*, . . . , i?Ar] and set J* = J*(p*) 
ENDFOR 



2N matrix summations are required. In summary, we need at most N+{R— 1)2N = N{2R— 1) 
matrix summations which is an upper bound for the complexity of greedy search. 



Algorithm 5 Greedy Bandwidth Distribution Algorithm 
Set Ro = = 0, . . . , i?7v = 0] and J° = Jf+p 
FOR i = l : R 
(A4.1) FOR k = l: N 

Add one bit to sensor k and compute det(J*(/c)) where 
= J^-i + AkiRk + 1) - AkiRk). 

ENDFOR 

(A4.2) Find the sensor p* for which det(J*(/c)) is the maximum: 

p* = arg max det{3^{k)). 

k 

(A4.3) Increment Rp* = Rp* + 1, update Rj = . . . , i?p*, . . . , Rn] and set J* = J*(p*). 
ENDFOR 
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IV. Simulation Results 

In this section, we illustrate the performance of different dynamic bandwidth distribution 
methods with numerical examples. For the convex relaxation problem, we solve the linear 
programming problem in (30) which is formulated to find the feasible initial point by using the 
"linprog" routine in MATLAB . The Newton method parameters e and r are selected according 
to [20]. Simulation results show that, for the convex relaxation scheme, the optimal solution is 
reached in around ten iterations. 

We evaluate the computation time of each bandwidth allocation approach by using the "etime" 
function of MATLAB averaged over 100 trials. In Fig. 3, the mean computation times of 
the considered suboptimal bandwidth allocation schemes are compared. Since the number of 
summations for A-DP increases linearly with N, for large number of sensors, the computation 
time of A-DP is less than the computation time of GBFOS and convex relaxation where the 
number of summations increase quadratically and cubically respectively. 





Convex Relax 




Approx. DP 




GBFOS 




Greedy 




20 30 40 50 

Number of Sensors 



Fig. 3: Computation time in seconds for convex relaxation, A-DP, GBFOS, and greedy search 
iR = 5). 



We next compare the MSB performances of the bandwidth distribution schemes based on 
optimal exhaustive search, A-DP, convex optimization, GBFOS and greedy search. In addition, 
we analyze the MSB performance of nearest neighbor bandwidth allocation, where all the 
bandwidth is assigned to the sensor which is nearest to the predicted target location. In our 
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simulations, we assume that sensors are grid deployed in a 6^ = 20 m x 20 m surveillance 
area as shown in Fig. 4-(a) and (b). We select Pq = 10^ and sensor observation noise cr^ = 1. 
The probability density function of the target's initial state, p(xo), is assumed to be Gaussian 
with mean /iq = [—8 — 8 2 2] and covariance Sq = diag[al 0.01 0.01] where we select 
ScTg = 2 so the initial point of the target remains in the ROI with very high probability. The target 
motion follows a white noise acceleration model and we consider two process noise parameters 
p = 2.5 X 10"'^ and p = 0.1. Measurements are assumed to be taken at regular intervals of 
V = 0.5 seconds and the observation length is 10 seconds. Namely, we perform target tracking 
over Ts = 20 time steps for each Monte-Carlo trial. The number of particles used in the particle 
filter is A'^, = 5000. We assume R = 5 bits of bandwidth is available at each time step for data 
transmission. The MSB at each time step is averaged over Tt^iais = 500 trials as. 



where in the v trial xj" and xj^ are the actual and estimated target states at time t respectively. 

In Fig. 4-(a) and (b), a WSN is illustrated where = 9 sensors track a target under the process 
noise parameters p = 2.5 x 10^'^ and p = 0.1 respectively. For p = 2.5 x 10^'^, the process 
noise is relatively small and the target trajectory is almost deterministic. For p = 0.1, the target 
trajectory has relatively large uncertainty. For the first time step of tracking. Table I presents 
each sensor's transmission probability for each quantization rate for the convex optimization 
based bandwidth allocation scheme with R = 5. Note that at t = 1, the target is relatively 
close to sensor 1 located at (—10 m.,— 10 m.). Then it is very likely that sensor 1 transmits 
its measurement using m = 5 bits because of the probability, gi 5 ^ 0.84. Rest of the sensors 
tend to remain silent since their transmission probabilities using bits are almost 1. As seen in 
Table II, the probabilistic transmission introduces a weak constraint on the bandwidth and on 
the average sensors transmit R bits to the fusion center. 

For = 9 sensors. Figs. 5-(a) and 5-(c) show the average number of sensors activated and 
Figs. 5-(b) and 5-(d) show the MSB at each time step of tracking averaged over 500 Monte- 
Carlo trials. Simulation results show that under p = 2.5 x 10"'^, convex relaxation, A-DP and 
GBFOS yield similar tracking performance to that of exhaustive search in terms of MSB. For 
p = 2.5 X 10""^, between the time steps 8 and 10, the target is relatively close to sensor 5 
located at (0,0). Hence, using exhaustive search, A-DP, convex relaxation, and GBFOS based 

October 26, 2011 DRAFT 




(41) 



JOURNAL OF MbX CLASS FILES, VOL. X, NO. X, SEPTEMBER 20XX 



25 



TABLE I: Transmission probabilities of each quantization rate for = 9 and i? = 5 at t = 1 
for the example illustrated in Fig. 4-(a). 





m = 


m — \ 


m = 2 


m — Z 


m = 4 


m = 5 


i = 1 


0.0021 


0.0010 


0.0011 


0.0037 


0.1482 


0.8440 


i = 2 


0.9877 


0.0082 


0.0003 


0.0017 


0.0012 


0.0009 


i = 3 


0.9895 


0.0053 


0.0021 


0.0014 


0.0010 


0.0008 


i = 4 


0.9906 


0.0031 


0.0031 


0.0017 


0.0011 


0.0003 


i = 5 


0.9888 


0.0057 


0.0023 


0.0014 


0.0010 


0.0008 


i = 6 


0.9895 


0.0053 


0.0021 


0.0014 


0.0010 


0.0008 


i = 7 


0.9895 


0.0053 


0.0021 


0.0014 


0.0010 


0.0008 


i = 8 


0.9895 


0.0053 


0.0021 


0.0014 


0.0010 


0.0008 


i = 9 


0.9895 


0.0053 


0.0021 


0.0014 


0.0010 


0.0008 



TABLE II: Mean and Standard deviation of the total number of transmitted bits using convex 
relaxation based bandwidth allocation method, R = 5. 





Mean 


Standard Deviation 


N = 9, p = 2.5x 10"^ 


5.0071 


1.9781 


AT = 25, p = 2.5 X 10"^ 


5.0291 


1.6264 


N = 9, p = 0.1 


5.0129 


1.3439 


N = 25, p = 0.1 


5.0093 


1.2119 



bandwidth allocation schemes, almost all the bandwidth is allocated to sensor 5. When the target 
is not relatively close to any of the sensors, as in time steps 2-6 and 12-17, the fusion center 
has relatively large uncertainty about the target location, so multiple sensors are activated with 
relatively coarse information which increases the estimation error as shown in Fig. 5-(b). After 
time step 17, the target approaches sensor 9 and by using exhaustive search, convex relaxation, 
A-DP and GBFOS, all the bandwidth is assigned to sensor 9 and then estimation error reduces 
again. The greedy bandwidth allocation scheme tends to activate more sensors all the time with 
relatively coarse information as compared to the other bandwidth allocation algorithms. With 
small process noise parameter (p = 2.5 x 10^'^), nearest neighbor based bandwidth allocation 
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becomes more accurate than greedy search since the target trajectory is highly deterministic and 
there is a small uncertainty on the predicted target location. However, the tracking performance 
of the nearest neighbor approach is still not as good as those for exhaustive search, convex 
relaxation, GBFOS, and A-DP. For p = 0.1, the uncertainty on target trajectory is relatively 
large and we observe a worse tracking performance as compared to the p = 2.5 x 10^^ case. On 
the other hand, still A-DP, convex optimization, and GBFOS perform equally well as exhaustive 
search in terms of MSE and outperform greedy search. For p = 0.1, nearest neighbor based 
bandwidth allocation introduces much larger estimation errors which are sometimes even greater 
than those obtained by the greedy search based dynamic bandwidth allocation scheme. 

For = 25 sensors. Figs. 6-(a) and 6-(c) show the average number of sensors activated and 
Figs. 6-(b) and 6-(d) show the MSE at each time step of tracking. Since the sensor density 
is increased, the bandwidth allocation schemes tend to assign all the available bandwidth to a 
single sensor which has more precise information about the target. This improves the tracking 
performance at each time step. For p = 2.5 x 10^'^ and p = 0.1 cases, convex relaxation, A-DP 
and GBFOS yield similar estimation performances and they significantly outperform the greedy 
search and nearest neighbor based bandwidth allocation approaches in terms of the MSE. 
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(b) 




(d) 



Fig. 5: N = 9, R = 5, Ttnai = 500 (a) Average number of active sensors, p = 2.5 x 10"^ , (b) 
MSB at each time step, p = 2.5 x 10^'^, (c) Average number of active sensors, p = 0.1, (d) MSB 
at each time step, p = 0.1. 



V. Conclusion 

In this paper, we studied the dynamic bandwidth allocation problem for target tracking in a 
WSN with quantized measurements. Under the bandwidth availability constraint, we proposed 
two bandwidth distribution schemes which are based on convex relaxation and approximate DP 
to maximize the determinant of the FIM. Simulation results show that convex relaxation, A- 
DP and GBFOS algorithms yield similar tracking performance, which is close to that provided 
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Time Step Time Step 

(c) (d) 

Fig. 6: = 25, i? = 5, Ttriai = 500 (a) Average number of active sensors, p = 2.5 x 10^^ , (b) 
MSB at each time step, p = 2.5 x 10""^, (c) Average number of active sensors, p = 0.1 (d) MSB 
at each time step, p = 0.1. 



by the optimal exhaustive search approach, and they outperform the greedy search and nearest 
neighbor approach significantly. Using the optimal solution of the convex optimization problem 
as the probability of transmission at each data rate, convex relaxation based bandwidth allocation 
satisfies the bandwidth constraint on the average while the other bandwidth distribution methods 
put a strict constraint on the bandwidth availability. In terms of computational complexity, A-DP 
is computationally more efficient than GBFOS and convex relaxation methods especially for a 
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large sensor network with a large A^. 

In this work, we developed and compared bandwidth allocation schemes in target tracking for 
one step ahead only. Our future work will cover extensions of proposed schemes to non-myopic 
scenarios. Multi-target tracking by dynamic bandwidth allocation will also be considered as a 
future research direction. 
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